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Compact objects, like neutron stars and white dwarfs, may accrete dark matter, and then be 
sensitive probes of its presence. These compact stars with a dark matter component can be modeled 
by a perfect fluid minimally coupled to a complex scalar field (representing a bosonic dark matter 
component), resulting in objects known as fermion-boson stars. We have performed the dynamical 
evolution of these stars in order to analyze their stability, and to study their spectrum of normal 
modes, which may reveal the amount of dark matter in the system. Their stability analysis shows 
a structure similar to that of an isolated (fermion or boson) star, with equilibrium configurations 
either laying on the stable or on the unstable branch. The analysis of the spectrum of normal modes 
indicates the presence of new oscillation modes in the fermionic part of the star, which result from 
the coupling to the bosonic component through the gravity. 



I. INTRODUCTION 

Scalar fields are of great interest in several fields of 
physics. In high energy physics they arise naturally in 
several unification theories, such as scalar-tensor theo- 
ries of gravitation from string theory. In cosmology, they 
have been considered to provide inflationary solutions in 
the early universe and an alternative explanation for dark 
energy. In addition, they have also been proposed as 
strong candidates of dark matter, the matter that is re- 
sponsible for the formation and evolution of structures 
in the Universe. For the latter type of models, one of 
the possibilities assumes that dark matter is composed 
by bosonic particles which may condensate into macro- 
scopic self- gravitating objects (i.e., self-gravitating Bosc- 
Einstcin condensates) commonly known as boson stars. 
Since the seminal paper in the late sixties by Ruffini and 
Bonazzola pj, boson stars in General Relativity have 
been extensively studied in may different contexts (for 
a recent review see @]). 

On the other hand, the formation of either a boson or 
a fermion star would be susceptible to subsequent mix- 
ture by fermions/bosons, and this fact opens a whole 
new possibility for the formation of objects made of both 
fermionic and bosonic particles. Even if one of these ob- 
jects is formed in a medium absent of either bosonic or 
fermionic particles, the latter may be accreted in later 
stages. In particular, bosonic dark matter particles may 
accrete on compact stars, and depending on the model 
considered, their effects on the star will be different and 
possibly observable. 

In the context of WIMPs, if the dark matter is self- 
annihilating, the released energy due to the WIMP an- 
nihilation inside the neutron star can increase the tem- 
perature and be observable in old stars Q . If it does not 
self-annihilate, the dark matter will cluster in a small re- 
gion at the center of the neutron star, increasing their 
compactness and ultimately leading to a gravitational 
collapse [J- Neutron stars may be therefore sensitive 
indirect probes of the presence of dark matter, and can 



be used to set constraints both on the density and on the 
physical properties of dark matter. 

Recent studies investigate possible changes in the 
structure of the star in the presence of dark matter, by 
using a two- fluid modelQ. In this paper, we perform 
a similar analysis by modeling systems which contain a 
fermionic compact star (we consider it to be a neutron 
star) and a bosonic dark matter component represented 
by a boson star. The resulting objects are known as 
fermion-boson stars. These mixed stars were first in- 
troduced by Henriques et. al. [1] (and further studied 
in @, Q ) , where the fermionic matter was described by a 
perfect fluid with the Chandrasekhar equation of state, 
while the bosonic component is modeled by using a real 
quantized scalar field as in Q. The bosons and the 
fermion particles are coupled only through gravity (no- 
tice however that non-minimal couplings with the scalar 
field can arise in other scenarios, such as in neutron stars 
with hidden extra dimensions [10( or in tensor-scalar the- 
ories of gravitation [Tl|). We will perform a dynamical 
analysis of these mixed stars by using a simple polytropic 
equation of state, as it is standard for cold neutron stars, 
and a complex scalar field to describe the bosonic com- 
ponent. 

The equilibrium configurations for cither an isolated 
boson or fermion stars are described, respectively, by the 
central value of the scalar field (f> c , and the central value of 
the fluid density p c 0, HH • These configuration are there- 
fore characterized by a single parameter o~ c , so in this case 
there are stability theorems [13l . Hij which indicate that 
the critical mass (separating the unstable from the stable 
branch) is located at the extrema dM/da c = 0. However, 
the mixed fermion-boson stars are parametrized not by 
one, but by two quantities (<p c , p c ). This implies that the 
analysis of stability is more complicated than in the iso- 
lated star case, since the previous stability theorems can 
not be directly applied. One can still analyze their stabil- 
ity, among other alternatives, by studying the radial per- 
turbations of these equilibrium configurations and then 
analyzing the eigenvalues of these modes in the linearized 
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equations as in |15l - ll9j , or by evolving dynamically these 
perturbations by solvingthe full non-linear equation of 
motion poU^ . In 13, 2f| Henriques et al. described a 
method to perform the analysis of stability of the boson- 
fermion stars by using the binding energy and the num- 
ber of bosonic and fermionic particles as a function of the 
two free parameters. In this paper, we propose a similar 
criterion, and our results are compared with the full nu- 
merical solution of the equations of motion. In addition 
to the stability analysis, we follow the migration of a star 
from an unstable to the stable branch, a process observed 
already in isolated boson and fermion stars. Finally, we 
study the dependence of the quasi-normal modes of the 
mixed star with respect to their fraction of bosonic mat- 
ter. 

The paper is organized as follows. In Sec. [II] we intro- 
duce the formalism used to obtain the set of evolution 
equations that describes the spacetime geometry and the 
boson-fermion matter contents. In Sec. IIIII we describe 
how to construct the initial data for the boson-fermion 
stars, and propose a method to find the stability of the- 
ses objects. The results of the dynamical evolution for 
equilibrium configurations (i.e., both stable and unsta- 
ble) are presented in Sec. lIVi together with the spectrum 
of the quasi-normal modes of the stable stars. Finally, 
conclusions and final remarks are presented in Sec. [V] 
Throughout this paper we use that the indices are a, b, . 
taken to run from to 3, while indices .. run from 
1 to 3. We also adopt the standard convention for the 
summation over repeated indices. 



II. FORMALISM 

Fermions minimally coupled to bosons can be modeled 
by considering a stress-energy tensor with contributions 
from a perfect fluid and a complex scalar field, in the 
form 



Tab 



y(fluid) 



T 



b ' 



(1) 



^r d) = [Po(l + e) + P}u a Ub + Pg ab , (2) 



[d a 4>* d b 4> + d a 4>d b 4>*} - 
g ab [d Q 0*d^ + m 2 |0p 



The perfect fluid is represented by the fermionic phys- 
ical (primitive) variables, namely the pressure P, rest- 
mass density p , internal energy e, and four- velocity u a , 
whereas the complex scalar field <j> describes a Bose- 
Einstcin condensate of bosonic particles of mass m. The 
fluid and the scalar field do not interact directly, and 
arc only coupled through gravity, as it is expected for 
WIMPS. The equations of motion for the fluid and the 
scalar field are obtained from the conservation laws of the 
stress-energy tensor and the baryonic number 



and the Klein-Gordon equation 
V a V> = m 2 



(4) 



(5) 



which together with the Einstein equations G a b = S^T a b 
constitute the system of equations governing the dynam- 
ics. 

We restrict our study to spherically symmetric stars, 
and then consider the time-dependent metric 

ds 2 = -a 2 (t,r) dt 2 + g rr (t, r)dr 2 + r 2 gee (t,r)dQ 2 . (6) 

The evolution equations for the spacetime are obtained 
by considering the Z3 formulation of the Einstein equa- 
tions [26j . which introduces the following independent 
quantities to form a first order system of equations, 



KJ 



d r a 
a 

1 d t g 



9' 



■d r g r 



-d r gee , 



2a g rr 



Ko 



1 dt. 



2a gee 



(7) 



The full system of equations for this formulation is in- 
cluded in appendix [X] The remaining freedom in the 
choice of coordinates of the line element $Q is related 
to the prescription for the lapse function, and a common 
option is the harmonic slicing condition 



dta = —a 2 trK . 



(8) 



where trK = K r r + 2Kg e . By using the metric ([5]), the 
equations of motion for the perfect fluid (jlj) and the scalar 



(3) field §5§ can be written explicitly as: 



dt(VlU) 



— d r (y/^av r D) ^/jav r D, 

-d r (VjaS r ) + ^a 
-d r (y/^jaS r r ) + ^ja 



S r r K r r + 2Se v K l 
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<t> t = (9 r (av^^r) + aVF ? [2(Pre + ^)0r + 2V^A^V t -m 2 5rr 



(9a) 
(9b) 
(9c) 
(9d) 
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where ^7 = ^/g rr ggg and we have introduced the auxil- 
iary fields 

4> r = d r (f>, <j) t = ^^-d t (f>, (10) 
a 

to reduce the Klein- Gordon equation to first order in 
space and time. 

The evolution of the fluid is described in terms of the 
conserved variables, namely the mass density D, the mo- 
mentum density S r and the energy density U. They are 
related to the primitive variables (i.e., the rest-mass den- 
sity p , the pressure P and the velocity v r ) by the fol- 
lowing relations 

D = poW, U = hW 2 -P, S r = hW 2 v r , (11) 

where h = po(l + e) + P is the enthalpy and W = 
— v r v r the Lorentz factor. In the right-hand-side of 
their evolution equations, the spatial projections of the 
stress-energy tensor take the form, 

S r r = hW 2 v r v r +P, S e e = P. 

During the evolution, the relations (jlip must be in- 
verted in order to obtain the primitive physical quantities 



(which are necessary for computing the rhs) from the con- 
served evolved fields. In general, this conversion can not 
be performed analytically, so appendix |B] explains in de- 
tail our numerical algorithm for obtaining the primitive 
fields. 



III. INITIAL DATA 

Initial data for the fermion-boson stars involves the in- 
trinsic metric gij and extrinsic curvature Kij on a given 
hyper-surface, as well as the fermionic fluid configura- 
tion in terms of its primitive variables (p, e,v l ) and the 
bosonic scalar field cf>. Assuming a static spherically sym- 
metric metric in Schwarzschild coordinates 

ds 2 = -a 2 {r)di 2 + a 2 {r)dr 2 + r 2 dtf , (12) 

a harmonic form of the scalar field (j>(t,r) = (j>(r)e~ lujt , 
and a star in hydrostatic equilibrium with v r = 0, the 
following system of ODEs is obtained: 



da 
dr 
da 
dr 
d4 
dr 

dr 
dP 
dr 



a f 1 



2 [r 
a { 1 



-(1 -a 2 ) +4nGr 



2 [r 

§(r), 



(a 2 - 1) +4irGr 



+ m 2 a 2 (j) 2 (r) + &(r) + 2a 2 p{\ + e) 



m 2 a 2 <p 2 {r) + $ 2 (r) +2a 2 P 



■\ if, 
r J a 2 (j> - [l + a 2 - ^Ga 2 r 2 (m 2 <l> 2 + p(l + e) - P)] — , 



-[p(l + e) + P]- 
a 



(13a) 
(13b) 
(13c) 
(13d) 
(13e) 



r 



The system is completed by choosing the equation of 
state (EoS) that relates the pressure with the other fluid 
quantities. As it is standard in simple models of cold 
stars, we will adopt here a polytropic equation of state 
P = Kp r , with the particular choice of T = 2 and 
K = 100, which corresponds to masses and compactness 
in the range of neutron stars (l2l |. 

We will use units such that c = 1, and the variables can 
be renormalizcd to absorb the factors G and to, so that 
the basic scale of the stars will be given by {K, V}. The 
final system is an eigenvalue problem for the frequency 
of the boson star uj as a function of two parameters; the 
central value of the scalar field (f> c and the density of the 
fluid p c . This system can be solved by using the Shooting 
Method HI. 



The appropriate boundary conditions for the scalar 
field and metric functions are obtained by imposing the 
conditions of regularity at the origin and asymptotic flat- 
ness at infinity. The condition at r = for the fluid pres- 
sure is obtained from the polytropic EoS as a function of 
p c . Thus, the full boundary conditions are 

o(0) = 1, o(0) = 1, <j>{0) = <pc, (14a) 
$(0) = 0, P(0) = Kp T c , (14b) 

(14c) 
(14d) 



lim a(r) = lim 



lim i 

r— too 



a(r) 

0, lim P(r) = 0. 



After the solution is found, a change of coordinates 
from Schwarzschild to maximal isotropic ones is per- 
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formed 



ds 2 



-a 2 (f)dt 2 + ^ 4 (f) (dr 2 + fW) , (15) 



which are more convenient for our numerical evolution 
and for future comparisons in three dimensions. All of 
our simulations will be shown in these coordinates, and 
for simplicity we will substitute f — > r hereafter. 

The total gravitational mass is computed by the 
asymptotic value of the metric coefficients 



r I 1 
M T = lim - ( 1 ff 



(16) 



The U(l) symmetry in the Lagrangian of the scalar field 
ensures the conservation of a Noether charge which can 
be associated with the number of bosons Nb [13, (HI- 
Correspondingly, the conservation of baryonic number al- 
lows to define a number of fermions Np. These quantities 
can be computed by integrating their densities, 



dN B 4TTau(j) 2 r 2 



dr 



dN, 



dr 



F = Anapr 2 



(17) 



Therefore, the radius of the fermionic/bosonic parts of 
the star can be defined as the surface containing 99% of 
the corresponding particles. 



A. Boson and fermion stars 

As a basic test of our initial data implementation we 
compare our equilibrium configurations with previously 
published results for isolated boson stars and fermion 
stars, which are the limits of our system of equations 
when p c — > and (p c — > respectively. 

Fig. ((T|) shows the total mass Mt of boson stars and 
fermion stars as a function of the corresponding radius 
Rgg. In agreement with the results of previous works, we 
have found that the maximum mass M max (i.e., the value 
of the mass that separates the stable M < M max from 
the unstable M > M ma x configurations) for the case of 
boson stars is M max = 0.633, whereas for fermionic stars 
is M max = 1.637 with r = 2 and K = 100. 



B. Mixted boson-fermion stars 

As we mentioned before, the equilibrium configurations 
of mixed boson-fermion stars are more involved and de- 
pend on the two parameters 4> c and p c - The total mass of 
the stars as a function of these parameters is plotted in 
Fig. [2] showing that the maximum mass is obtained for 
the isolated neutron star case (i.e., when cf> c = 0). This is 
direct consequence of our choice of the parameters {K, T} 
in the equation of state, that sets the scale and the com- 
pactness of the mixed stars. With the current choice, 
the stars will be composed predominantly by fermions, 
which can produce stars with much higher compactness 
than boson stars. 




FIG. 1: Initial data of isolated stars The total masses of the 
boson Mb and fermion Mf stars, as functions of their cor- 
responding radius The maximum mass agrees in each 
case with previous results found in the literature, namely 
M max = 0.633 for boson stars, and M max = 1.637 for fermion 
stars. 
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FIG. 2: Initial data of mixed fermion-boson stars The total 
mass of the equilibrium configurations of the mixed stars, as 
a function of <f> c and p c - The maximum mass for a given value 
of p c , is always found when <j} c = 0. 



The profiles of the different non-trivial fields for a rep- 
resentative case are plotted in Fig. [3] , which clearly sat- 
isfy the regularity conditions at the origin and asymptotic 
flatness. The presence of the fermionic fluid produces 
a deeper gravitational potential than the one produced 
solely by the boson star, therefore contracting the bosonic 
component to a smaller radius, comparable to the one of 
the fermionic matter. 

For a fixed value of the total mass Mr, we find that the 
number of bosons increases for (f> c > 0, reaches a max- 
imum, and then decreases. Notice that since the mass 
is kept fixed, the central density p c must change as we 
vary <p c . The number of fermions has consequently the 
complementary behavior; it decreases until reaching a 
minimum and then increases. The same profiles are ob- 
served in these quantities, switching Np by Nb, when 
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FIG. 3: Initial data of mixed fermion-boson stars. The pro- 
files of the scalar field <j>{r), the fermionic density p(r), and 
the conformal factor ip(r) for one typical configuration corre- 
sponding to Nb/N f = 0.1 and M T = 1.4. 



they are represented MS c\ function of p c instead of </> c . 
This behavior is illustrated in Fig. |4l where the num- 
ber of particles is plotted as a function of <j) c and p c for 
the configurations with mass Mr = 1.4. In the next sec- 
tion, we will describe the evolution of the two equilibrium 
configurations marked with symbols in Fig. [4l one with 
Nb = !0%Np which is on the left of the maximum value, 
and the other with Nb = 10.7%Np which is on its right. 
Our simulations show that the first configuration is in the 
stable branch, while the second one lies in the unstable 
branch, indicating that the maximum/minimum of these 
curves is the critical point. 

The stability analysis for the boson-fermion stars, 
which can be performed for instance by solving the per- 
turbed equations, is much more complicated than for iso- 
lated boson or fermion stars. The main reason, as it was 
mentioned before, is that these mixed configuration have 
two free parameters (i.e., the central values of the scalar 
field 4> c and the density of the perfect fluid p c ) instead of 
just one, so that the stability theorems for a single pa- 
rameter solutions can not be directly applied. One possi- 
bility would be to find a method to perform the stability 
analysis such as those described in [3, H, [HI, HH • 
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FIG. 4: Initial data of mixed fermion-boson stars The number 
of fermions Nf and bosons Nb for the equilibrium configura- 
tions as a function of (f> c (top panel) and p c (bottom panel). 
The position of the maximum/minimum corresponds to the 
critical point which separates the stable and the unstable so- 
lutions. The two configurations considered in the next section 
are marked, one on each side of the maximum/minimum, cor- 
responding to N B = 10%N F (t) and Nb = 10.7%N F (If). 



Instead of performing this stability analysis, we pro- 
pose a different way to define the stability of the mixed 
configurations based on the critical value of the num- 
ber of particles. Our criterion states that the configura- 
tions with the number of bosons (fermions) on the left 
of the maximum (minimum) are stable configurations, 
while configurations that are on the right of the maxi- 
mum (minimum) are unstable. The solution space of sta- 
ble/unstable configurations using this criterion is shown 
in Fig. [5] This will be validated through our numerical 
simulations in the next section. 
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FIG. 5: Initial data of mixed fermion-boson stars. Regions 
of stability /instability for the equilibrium configurations of 
the mixed boson-fermion stars, according to the criterion of 
maximum/minimum of the number of bosons/fermions for a 
fixed Mr, see the text for more details. We mark the two 
configurations corresponding to Nb = 10%JVi?(|) and Nb = 
10.7%N F (1). 



IV. NUMERICAL SIMULATIONS 

In this section we analyze the dynamics of mixed stars, 
and address different issues like the stability of these sys- 
tems or their spectrum of normal modes. In order to 
determine the properties of the mixed star equilibrium 
configurations described in the previous section, we per- 
formed long-term numerical evolutions of the discretizcd 
Einstein-Klcin-Gordon-Hydrodynamic system ([9]). 

First of all, we write the system in flux conservative 
form 



d t \5 + d k F k {V) = S(U). 



(18) 



so that we can apply numerical algorithms based on 
Finite Volume methods. The spatial discretization of 
the geometry and the boson fields is performed using a 
third order accurate Finite Volume method [2f|, which 
can be viewed as a fourth order finite difference scheme 
plus third order adaptive dissipation. The dissipation 
coefficient is given by the maximum propagation speed 
in each grid point. For the fluid matter fields, we 
use a High Resolution Shock Capturing method with 
Monotonic-Centered limiter. The time evolution is per- 
formed through the method of lines using a third order 
accurate Strong Stability Preserving Runge-Kutta inte- 
gration scheme [28[, with a Courant factor of At/Ar = 
0.25 so that the Courant-Fricdrichs-Lcvy (CFL) condi- 
tion dictated by the principal part of the equations is 
satisfied. Most of the simulations presented in this work 
have been done with a spatial resolution of Ar = 0.01, 
in a domain with outer boundary situated at r = 600. 
We use maximally dissipativc boundary conditions for 
the spacetime variables and the boson fields, and outflow 
boundaries for the fluid matter fields. 
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FIG. 6: Evolution of stable fermion-boson stars The central 
values of the density and the (peaks of the oscillatory) scalar 
field (top), and the mass and the number of bosonic and 
fcrmionic particle (bottom). All the quantities remain very 
close to their initial values, suggesting that the star is stable 
against perturbations. 



We start with the dynamical evolution of the mixed 
equilibrium solution corresponding to Nb = 10%Nf and 
total mass Mt = 1-4, shown in Fig. 2) Since it is located 
on the left of the critical value, we expect this configura- 
tion to be stable. 

The evolution displays a combination of the behav- 
iors that are typical for isolated boson and fermion stars. 
The scalar field oscillates with its characteristic eigenfrc- 
quency, while the fluid density oscillates slightly around 
its initial state due to the perturbation introduced by 
the numerical truncation errors. The values of the peaks 
of the oscillatory scalar field 0™ aa: and the fluid density 
po at the center of the star are plotted as a function of 
time in the top panel of Fig. |6l while the total mass Mt 
and the number of particles Nb , Np are displayed in the 
bottom panel. 

These quantities remain very close to their initial value 
for many dynamical times (except for a tiny drift due to 
numerical dissipation), indicating that the configuration 
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is indeed stable. In order to asses the robustness and ac- 
curacy of our numerical implementation, we have evolved 
this configuration with three different spatial resolutions 
Ar = (0.02,0.01,0.005), in a domain of r = 600, for 
t ps 2000, finding that the numerical solution converges 
at second order. The energy constraint (|A6j) is small dur- 
ing the evolution and converges to zero, as it is shown in 

Fig.m 
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FIG. 7: The energy constraint (|A6[) for three different reso- 
lutions Ar = 0.005 (red), Ar = 0.01 (green), and Ar = 0.02 
(blue), showing second order convergence. 
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B. Unstable stars 



FIG. 8: Evolution of unstable fermion-boson stars Same as 
fig. [6] for a star on the right of the critical curve. The central 
values of the density and (the peaks of) the scalar field depart 
quickly from their initial values, indicating that the star is 
unstable. The evolution becomes non-linear and describes the 
migration of the star from the unstable to the stable branch. 



The numerical evolution of the equilibrium configura- 
tion with Nb = 10.7%-ZVf and Mt — 1-4 presents a more 
dynamical behavior. This configuration lies on the right 
of the critical values of the number of particles Np and 
Nb in Fig. 01 indicating that it must be unstable under 
perturbations. 

The initial stage of the evolution is similar to the pre- 
vious case of a stable star, with the scalar field oscillat- 
ing mainly with its eigenfrequency, and the neutron star 
oscillating due to the perturbation introduced by the in- 
herent numerical truncation errors. However, these oscil- 
lations grow rapidly in amplitude, driving the dynamics 
to a non-linear regime; the star is eventually migrating 
from the unstable to the stable branch. The central val- 
ues po and the maximum </>™ aa: , the total mass Mt, and 
the number of particles Nb and Np, are plotted in Fig. [5] 
The central values show large variations until t ps 10000, 
then change slowly, and finally settle down onto a new 
stable configuration with practically the same number of 
bosonic particles, but with around 10% less mass and 
number of fermions than the original star. 



C. Quasi-Normal Modes of the stable stars 

As it has already been mentioned, the fermion-boson 
star will oscillate around its stable configuration due to 
the perturbations introduced by the numerical trunca- 
tion errors, in a similar way as an isolated fcrmion or 
boson star. These perturbations will excite the charac- 
teristic modes of the mixed star, so that the oscillations 
will be a superposition of normal modes, each one with 
a characteristic frequency. 

By analyzing the central oscillations of the different 
fields, and in particular, of the central density of the 
star, we can study the structure of the normal modes of 
the fermion-boson stars. The frequencies of the normal 
modes are well-known for both isolated neutron and bo- 
son stars, but they have not been yet studied for mixed 
stars. The pulsations of compact objects are of great 
importance for relativistic astrophysics, because they of- 
fer the possibility of extracting information about the 
star (for instance the radius, mass and equation of state) 
from the detection of the associated gravitational waves 
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(see [29[ for a review). Although our spherical symmetry 
assumption only allow us to study radial modes (i.e., the 
fundamental mode and its overtones), it is still represen- 
tative to show how these modes may change in a neutron 
star in the presence of a bosonic dark matter component 
that couples to fermions only through gravity. 



We will restrict our analysis to a fermion-boson star 
with total mass Mj = 1.4, and parametrize different 
mixed stars by increasing the amount of bosons relative 
to fermions, corresponding to the fractions Ns/Np = 
{0, 2.5, 5, 7.5, 10}%. The details of the parameters of the 
stars are summarized in TablclU We have evolved for long 
times t « 6000 in order to get at least 50 oscillations of 
the central density, which will produce a clear spectrum 
with sharp peaks in the frequency domain. The Fourier 
transform of this quantity is shown in the top panel of 
Fig. IH1 As an additional check of our code, we can com- 
pare the known frequencies of the fundamental mode and 
its overtones for a (fermion-only) neutron star (as com- 
puted either by using perturbation theory or numerical 
evolutions, see for instance [30l ]) with the ones obtained 
from our simulation for the purely fermionic case (corre- 
sponding to the circles on the left in the bottom panel 
of Fig. [5]). The difference is always smaller than 1%, 
confirming the accuracy and correctness of our results. 



We now turn our attention to the boson-fermion case. 
The fundamental mode, which is usually a function of 
the mean density of the star, remains roughly constant 
except for the largest boson fraction, for which it shifts 
slightly towards smaller frequencies. The overtones, at 
higher frequencies, display more interesting features. In 
particular, there appears to be a splitting of these modes. 
The original neutron star overtones, displayed with cir- 
cles in Fig. [SJ are the dominant ones for small number of 
bosons. The power of the new oscillation modes increases 
with the boson fraction, suggesting that their origin is 
the gravitational coupling with the scalar field. The fre- 
quency of the overtones has a significant drift towards 
higher values as the fraction of bosons increases. 



The main features of this spectrum can be qualitatively 
explained in a very simple way. The new quasi-normal 
modes, which were not present for isolated fermionic 
stars, corresponds to the quasi-normal modes of the bo- 
son star. The oscillations in the bosonic part propagate 
to the fermions through gravity. As the fraction of bosons 
increases, so does the relative importance of the scalar 
field with respect to the fluid density, producing the ob- 
served growth in the amplitude of these modes. Conse- 
quently, the splitting effect is just a superposition of the 
quasi-normal modes of the boson and the neutron star. 
The drift in the frequencies is an effect of the change in 
radius and mean-density of the star as the fraction of 
boson changes. 
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FIG. 9: Normal modes of the fermion-boson stars (Top) 
Fourier spectrum of the central value of the fluid den- 
sity po for several stable configuration with Nb/Nf = 
{0, 2.5, 5, 7.5, 10}% and M T = 1.4. (Bottom) Frequencies cor- 
responding to the first, second and third modes of the isolated 
neutron star, as a function of the boson fraction. Notice the 
appearance of new oscillation modes, not present for an iso- 
lated neutron star. See the text for more details. 



V. CONCLUDING REMARKS 

We have studied in some detail the numerical evolu- 
tion of equilibrium configurations of mixed boson-fermion 
stars. Our results confirm the existence of stable and un- 
stable branches of equilibrium configurations. We also 
defined a stability criterion based on the variation of the 
number of bosonic and fermionic particles, for a given 
fixed value of the total mass, as a function of the central 
values of the scalar field amplitude and the fluid density. 
This criterion states that the equilibrium configurations 
located on the left of the maximum (minimum) number 
of bosons (fermions) are stable, whereas the configura- 
tions located on the right are unstable. We were able 
to determine the curve that separates the stable branch 
from the unstable one, in the plane formed by the cen- 
tral values of the scalar field cj) c and fluid density p c . We 
also verified that the correct solutions are obtained in the 
limiting cases of an isolated boson or fermionic star, by 
comparing with the results of previous studies. 

In order to assess the stability criterion, we performed 
the numerical evolution of the fully non-linear equations 
of motion for two types of solutions. For the stable con- 
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TABLE I: Properties of the fermion-boson star models used in the simulations. All the stars have a total mass Mt — 1.4. 
The columns report: the fraction of boson particles, the central value of the scalar field <f> c , the central density p c , the internal 
frequency of the scalar field lob, the total radius of the star Rt, the number of bosonic particles Nb, and the radius of the 
bosonic and the fermionic components, Rb and Rf, respectively. Notice that the largest fraction of Nb/Nf for a stable 
configuration is reached for the maximum value of Nb and the minimum value of Nf, but the precise value of the fraction 
depends on the value of the total mass (in the present case we get Nb/Nf ~ 13%); larger ratios Nb/Nf can be obtained for 
smaller values of Mt ■ 



figuration, the central values of the scalar field and the 
fermionic density remain constant in time during the nu- 
merical evolution, while the unstable star migrates to a 
stable configuration by ejecting out some of the initial 
mass. 

We also studied the structure of the normal modes and 
overtones of these mixed stars by performing long term 
numerical evolutions for configurations with a fixed to- 
tal mass but with different boson to fermion ratios. As 
expected, new oscillation modes appear in the frequency 
spectrum of the stars, when compared to the fermion- 
only case; the appearance of the new overtones is justi- 
fied because of the gravitational coupling of the fermionic 
perfect fluid with the scalar field, which has its own os- 
cillation modes. 

As we mentioned before, an accurate classification of 
the properties of boson- fermion stars is necessary in order 
to investigate the possible existence of bosons trapped in- 
side, for instance, in neutron stars. One possible indica- 
tion of such a phenomenon would be the frequency shifts 
and splitting in the vibration spectrum of the stars. 

Another case of astrophysical interest is the possible 
existence of bosonic dark matter in galactic halos, an 
idea that has drawn some attention in the specialized 
literature in recent years 0, [3ll - l4l| . This type of mod- 
els refers to the other extreme case, in which the star 
is dominated by the bosonic component. In the context 
of galaxies, these mixed fermion-boson stars could model 
the galaxy halo with a boson star, and the gas with a 
fermionic component. Our present results indicate that a 
boson-dominated galaxy halo must keep its stability fea- 
tures after the inclusion of fermions; however, more work 
is needed to determine the properties of the equilibrium 
configuration that may be detected through astrophysi- 
cal observations. This is work in progress that we expect 



to report elsewhere. 
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Appendix A: The evolution system of equations 

We consider the Z3 formulation of the Einstein equa- 
tions in spherical symmetry [23j as the evolution system 
for the space-time geometry. The system is regularized 
at the origin using the following transformation of the 
momentum constraint: 



which ensures the cross-cancellation of the factors 1/r 
in the fluxes, and 1/r 2 in the sources. The sources still 
have terms like 1/r times other variables that contain ra- 
dial derivatives of the metric coefficients. However, these 
terms do not create problems at r — > 0, as the radial 
derivatives of any diffcrentiablc function must vanish at 
the origin. Thus, the equations of motion read 
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d t A r = -d/atrK], 
dtD rr r = — d r [aK r r ] , 
d t D re 9 = -d r [aK e e ], 

8 t Z r = -d r [2aK e e ] + 2a[{K/ - Kg e )(p r0 e + -) - K r ' 

+A r K 6 e + ^{Kg e - K/) - AnSr} , 
Ar Qgg J 



d t K 9 6 



At 



d t K/ = -d r 



Ar gee 
2 

ag rr (A r + -D r 



\z r )] + a{{K/) 2 + ^Kg 6 {K r 



1 - 



Kg 



-g rr D r /A r + -[g rr {D r / - A r - AZ r ) + g ee (D r g e - A r )} 
Sr 

9ee' 



\g rr (D r g° + 1) {D r / A r ) + SttQ - + S °) } , 

-a{iAV(-AV r +4^ e e ) 



-d r 



ag 



1 fl 2 
3^+3*' 



--[<T(A r - 2£> r / - 4Z r ) + .g e0 (A r - 2£> r fl fl )] 
or 



"3 S 



Ar V Qgg 



(D rr r ~ D r g" - 2A r ) 
\g" (Dj + \) {Dr/ - AA r ) + 8tt(1 - ^ + S e e ) } , 



(Ala) 
(Alb) 
(Ale) 



(Aid) 



(Ale) 



(Alf) 



where Z r is the vector associated with the Z3 formula- 
tion, and trK = K r r + 2K e g is the trace of the extrinsic 
curvature. In Sec. |IT]we defined the matter terms of the 
fcrmionic fluid {D, U, S r , S/ , Sg B }, and the auxiliary 
variables {A r , D r / , D r g 6 } which we introduced in order 
to reduce the full system in Eqs. HI and IA11 to first 
order in space and time. 

The total matter terms are given by 

r = ^(g^^ + g^^r + Vm + U, (A2a) 

S r = — [y/S^^ + V^M^ + Sr, (A2b) 

S/ = ±[g rr r t <t>t + g rr ct>*r<t>r-V(ct>)] + S/,(A2c) 

Se e = l -[g rr ^t-g rr rAr~Vm+~Sg\{A2d) 

and the total mass of the mixed stars is calculated from 
the Tolman mass defined as 

M T = j \t ° ~TC)^dx\ (A3) 

= Att J r 2 ayfg~/ggg{T + S/ + 2S e e )dr . 

On the other hand, the number of fcrmionic particles 
associated to the mass of the fermionic fluid given by 



Np = Att r z a^g—.ggg{U + S/ + 2Sg u )dr . (A4) 



The number of bosonic particles can be associated to 
the Noether charge [l[ of the scalar field, which can be 
computed as 



N 



B 



An / r 2 y/gr~r~geo((l>*dt<t> - <f>dt<j>*) / y/cfrrdr . (A5) 



The Hamiltonian constraint takes the form 

2 f „ ^ ^ a „ ^ f) ( ^ a 2 

g rr 



H = —\-2diD r g 9 -3D r g e (D rt 

Qrr y V r 

(1 - grrg 6e ) 



+g rr K 9 "{Kg« + 2K/) 



-2D r / \ - + D r 



%TVg rr T^. 



(A6) 



Appendix B: The transformation from conserved to 
primitive quantities 

From the definition of the conserved quantities 

D = p W, U = hW 2 ~P, S r = hW 2 v r , (Bl) 

one obtains the primitives {p,P,v r , e} after each time 
integration of the equations of motion. This is not trivial, 
mainly because the enthalpy h = p(l + e) + P, and the 
Lorentz factor W = 1/\A — v r v r , are defined as func- 
tions of the primitives. 

We are adopting a recovery procedure which consists 
in the following steps: 
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1. From the first thermodynamics law for adiabatic 
processes, it follows that 

P = (T-l)pe. (B2) 

Substituting the definition of the entalphy in the 
equation of state above, we write the pressure as 
a function of the conserved quantities and the un- 
known variable x — hW 2 . 

2. Using the previous step, the definition of U be- 
comes: 

U = hW 2 - P 

= hW 2 - { -^{h-p) 

= W2(l-£_l)+£_l p , (B3) 

where T is the adiabatic index corresponding to an 
ideal gas. 

3. Then, the function 



must vanish for the physical solutions. The roots 
of the function f(x) = can be found numerically 
by means of an iterative Newton-Raphson solver, 
so that the solution at the n + 1-iteration can be 
computed as 

f( x n) /rir\ 

X n +1 = X n - — r, (B5) 

J \ x n) 

where f'(x n ) is the derivative of the function f(x n ). 
The initial guess for the unknown x is given in the 
previous time step. 

4. After each step of the Newton-Raphson solver, we 
update the values of the fluid primitives as 

p = D/W, P = x-U, v r = S r /x, (B6) 

where W 2 = x 2 /(x 2 - S r S r ). 

5. Iterate steps 3 and 4 until the difference between 
two successive values of x falls below a given thresh- 
old value of the order of 1CU 10 . 
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